/*
* INSANE - Interactive Structural Analysis Environment
*
* Copyright (C) 2003-2005
* Universidade Federal de Minas Gerais
* Escola de Engenharia
* Departamento de Engenharia de Estruturas
* 
* Author's email :    insane@dees.ufmg.br
* Author's website :  http://www.dees.ufmg.br/insane
* 
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or any later version.
* 
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
* 
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/
package br.ufmg.dees.insane.shape;

import br.ufmg.dees.insane.util.IMatrix;
import br.ufmg.dees.insane.util.IVector;

/** A class representing shape referring to the
 *  rectangular plate-bending element nonconforming of 
 *  four node developed by Melosh, Zienkiewicz and Cheung,
 *  called MZC.
 * @author Samir
 * @since 01/06/2006
 */
public class MZC extends ShapeC1 {

	private static final long serialVersionUID = 1L;

	/* (non-Javadoc)
     * @see br.ufmg.dees.insane.shape.Shape#getShapeFunction(double[], IMatrix)
     */
	public IVector getShapeFunction(double[] g, IMatrix cN)
    {
		IMatrix J = this.getJacobian(getDerivedShapeFunction(g,cN),cN);
		double a = this.getCte(Math.abs(J.getElement(0,0)),Math.abs(J.getElement(1,0)));
		double b = this.getCte(Math.abs(J.getElement(0,1)),Math.abs(J.getElement(1,1)));
		
        double x = g[0]; // x = xi = (x-xo)/a
        double e = g[1]; // e = eta = (y-yo)/b
        
        IVector n = new IVector(12);
        
        double n0 = 0.125*((1-x)*(1-e)*(2-x-e-(x*x)-(e*e)));
        double n1 = 0.125*((x*x)-1)*(x-1)*(1-e)*a;
        double n2 = 0.125*((e*e)-1)*(e-1)*(1-x)*b;

        double n3 = 0.125*((1+x)*(1-e)*(2+x-e-(x*x)-(e*e)));
        double n4 = 0.125*((x*x)-1)*(x+1)*(1-e)*a;
        double n5 = 0.125*((e*e)-1)*(e-1)*(1+x)*b;

        double n6 = 0.125*((1+x)*(1+e)*(2+x+e-(x*x)-(e*e)));
        double n7 = 0.125*((x*x)-1)*(x+1)*(1+e)*a;
        double n8 = 0.125*((e*e)-1)*(e+1)*(1+x)*b;

        double n9 = 0.125*((1-x)*(1+e)*(2-x+e-(x*x)-(e*e)));
        double n10 = 0.125*((x*x)-1)*(x-1)*(1+e)*a;
        double n11 = 0.125*((e*e)-1)*(e+1)*(1-x)*b;
        
        n.setElement(0,n0);
        n.setElement(1,n1);
        n.setElement(2,n2);
        n.setElement(3,n3);
        n.setElement(4,n4);
        n.setElement(5,n5);
        n.setElement(6,n6);
        n.setElement(7,n7);
        n.setElement(8,n8);
        n.setElement(9,n9);
        n.setElement(10,n10);
        n.setElement(11,n11);
        
        return (n);
    }

    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.shape.Shape#getDerivedShapeFunction(double[], IMatrix cN)
     */
	/* Utilizado as 1as. derivadas do elemento Q4, uma vez que este elemento
	 * deve ser utilizado somente na forma retangular.
	 */
	public IMatrix getDerivedShapeFunction(double[] g, IMatrix cN) {
		double a = g[0];
		double b = g[1];
		
		IMatrix dl = new IMatrix(2,4);
		
	//derivadas em relacao a xi (a = xi e b = eta)
		double dl0 = ((0.25)*(-1+b));
		double dl1 = ((0.25)*(1-b));
		double dl2 = ((0.25)*(1+b));
		double dl3 = ((0.25)*(-1-b));
		
	//derivadas em relacao a eta (a = xi e b = eta)
		double dl9 = ((0.25)*(-1+a));
		double dl10 = ((0.25)*(-1-a));
		double dl11 = ((0.25)*(1+a));
		double dl12 = ((0.25)*(1-a));
		
	//derivadas em relacao a xi na linha 0 da IMatrix
		dl.setElement(0,0,dl0);
		dl.setElement(0,1,dl1);
		dl.setElement(0,2,dl2);
		dl.setElement(0,3,dl3);
		
	//derivadas em relacao a eta na linha 1 da IMatrix
		dl.setElement(1,0,dl9);
		dl.setElement(1,1,dl10);
		dl.setElement(1,2,dl11);
		dl.setElement(1,3,dl12);
		
		return (dl);
	}

	/* (non-Javadoc)
     * @see br.ufmg.dees.insane.shape.Shape#getSecondDerivedShapeFunction(double[], IMatrix)
     */
	public IMatrix getSecondDerivedShapeFunction(double[] g, IMatrix cN) {

		IMatrix J = this.getJacobian(getDerivedShapeFunction(g,cN),cN);
		double a = this.getCte(Math.abs(J.getElement(0,0)),Math.abs(J.getElement(1,0)));
		double b = this.getCte(Math.abs(J.getElement(0,1)),Math.abs(J.getElement(1,1)));
		
		double x = g[0]; // x = xi = (x-xo)/a
        double e = g[1]; // e = eta = (y-yo)/b
        
        IMatrix dll = new IMatrix(3,12);
        
        double dll0 = 0.75*x*(e-1)/(a*a);
        double dll1 = 0.75*e*(x-1)/(b*b);
        double dll2 = 0.125*(-4+3*(x*x)+3*(e*e))/(a*b);
        
        double dll3 = 0.25*(-3*x+3*(x*e)+1-e)/a;
        double dll4 = 0;
        double dll5 = 0.125*(3*(x*x)-2*x-1)/b;
        
        double dll6 = 0;
        double dll7 = 0.25*(-3*e+3*(x*e)+1-x)/b;
        double dll8 = 0.125*(3*(e*e)-2*e-1)/a;
        
        double dll9 = -0.75*x*(e-1)/(a*a);
        double dll10 = -0.75*e*(1+x)/(b*b);
        double dll11 = -0.125*(-4+3*(x*x)+3*(e*e))/(a*b);
        
        double dll12 = 0.25*(-1+e-3*x+3*x*e)/a;
        double dll13 = 0;
        double dll14 = 0.125*(2*x+3*(x*x)-1)/b;
        
        double dll15 = 0;
        double dll16 = -0.25*(3*e+3*(x*e)-1-x)/b;
        double dll17 = -0.125*(3*(e*e)-2*e-1)/a;
        
        double dll18 = 0.75*x*(1+e)/(a*a);
        double dll19 = 0.75*e*(1+x)/(b*b);
        double dll20 = 0.125*(-4+3*(x*x)+3*(e*e))/(a*b);
        
        double dll21 = -0.25*(1+e+3*x+3*(x*e))/a;
        double dll22 = 0;
        double dll23 = -0.125*(2*x+3*(x*x)-1)/b;
        
        double dll24 = 0;
        double dll25 = -0.25*(1+x+3*e+3*(x*e))/b;
        double dll26 = -0.125*(2*e+3*(e*e)-1)/a;
        
        double dll27 = -0.75*x*(1+e)/(a*a);
        double dll28 = -0.75*e*(x-1)/(b*b);
        double dll29 = -0.125*(-4+3*(x*x)+3*(e*e))/(a*b);
        
        double dll30 = -0.25*(3*x+3*(x*e)-1-e)/a;
        double dll31 = 0;
        double dll32 = -0.125*(3*(x*x)-2*x-1)/b;
        
        double dll33 = 0;
        double dll34 = 0.25*(-1+x-3*e+3*(x*e))/b;
        double dll35 = 0.125*(2*e+3*(e*e)-1)/a;
        
    	dll.setElement(0,0,dll0);	dll.setElement(0,1,dll3);	dll.setElement(0,2,dll6);	dll.setElement(0,3,dll9);
    	dll.setElement(1,0,dll1);	dll.setElement(1,1,dll4);	dll.setElement(1,2,dll7);	dll.setElement(1,3,dll10);
    	dll.setElement(2,0,dll2);	dll.setElement(2,1,dll5);	dll.setElement(2,2,dll8);	dll.setElement(2,3,dll11);

    	dll.setElement(0,4,dll12);	dll.setElement(0,5,dll15);	dll.setElement(0,6,dll18);	dll.setElement(0,7,dll21);
    	dll.setElement(1,4,dll13);	dll.setElement(1,5,dll16);	dll.setElement(1,6,dll19);	dll.setElement(1,7,dll22);
    	dll.setElement(2,4,dll14);	dll.setElement(2,5,dll17);	dll.setElement(2,6,dll20);	dll.setElement(2,7,dll23);

    	dll.setElement(0,8,dll24);	dll.setElement(0,9,dll27);	dll.setElement(0,10,dll30);	dll.setElement(0,11,dll33);
    	dll.setElement(1,8,dll25);	dll.setElement(1,9,dll28);	dll.setElement(1,10,dll31);	dll.setElement(1,11,dll34);
    	dll.setElement(2,8,dll26);	dll.setElement(2,9,dll29);	dll.setElement(2,10,dll32);	dll.setElement(2,11,dll35);
        
        return (dll);
    }

    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.shape.Shape#getNaturalNodalCoords()
     */
    public IMatrix getNaturalNodalCoords() {
    	IMatrix natCoords = new IMatrix(2,4);
    	
    	//coordenadas naturais de xi
    		natCoords.setElement(0,0,-1);
    		natCoords.setElement(0,1,1);
    		natCoords.setElement(0,2,1);
    		natCoords.setElement(0,3,-1);
    		
    	//coordenadas naturais de eta
    		natCoords.setElement(1,0,-1);
    		natCoords.setElement(1,1,-1);
    		natCoords.setElement(1,2,1);
    		natCoords.setElement(1,3,1);
    		
        
        return natCoords;
    }

    /* (non-Javadoc)
     * @see br.ufmg.dees.insane.shape.Shape#getType()
     */
    public String getType() {
        return "MZC";
    }
    
    /** Return the jacobian matrix.
	*@param dl The matrix of local derivates.
	*@param cN The matrix of nodal coordinates.
	*@return dl * cN The jacobian matrix.
	*/
	private IMatrix getJacobian(IMatrix dl, IMatrix cN)
	{
		IMatrix jacobian = new IMatrix(dl.getNumRow(),cN.getNumCol());
		jacobian.setZero();
		jacobian.mul(dl, cN);
		return(jacobian);
	}
	
	private double getCte(double x, double y){
		if(x > y){
			return x;
		}else{
			return y;
		}
	}
}

